本学期总共六个实验,先后完成了对酵母YJM320的全基因组的模拟测序、序列组装、同源搜索、从头预测、启动子(TATA box)预测、全基因组数据可视化。
摘要
本学期总共六个实验,先后完成了对酵母YJM320的全基因组的模拟测序、序列组装、同源搜索、从头预测、启动子(TATA box)预测、全基因组数据可视化。
具体内容如下:利用art illumina工具对原始的酵母YJM320的全基因组进行模拟建模,利用fastqc对建模结果进行质控分析(质控结果很好)。然后对建模结果利用SOAPdenovo进行序列组装,组装结果为Scaffold265个,其N50 258107bp;contig大于100bp的共43082个,其N50为273bp。结果利用blastn和quast同原始基因组进行比对,两种方法计算覆盖度分别为75.94%和87.6%,利用IGV工具观察基因组,发现在12号染色体中有一段rRNA重复序列无法组装,这也是影响最终覆盖率的主要原因。
分别采用同源搜索(利用tblastn比对)和从头测序(Augustus)的方法构建全基因组,对所得结果同原基因组进行gffcompare比对,结果为从头测序结果优于同源建模。IGV是数据可视化的工具,利用它可以查看同源搜索和从头预测的结果。
采用HMM结合bootstrap预测全基因组中的TATA box分布,绘制打分图和ROC曲线图,有70%匹配的TATA box下游5000bp内有发现基因。
一、材料和方法
1.1、分析平台
1.1.1 硬件平台
(1).CPU: Intel(R) Core(TM) i7-7700HQ CPU @ 2.80GHz
(2). 内存:8 GB 1600 MHz DDR3
(3). 硬盘:磁盘 0 (C:) TOSHIBA THNSNK128GVN8 M.2 2280 128GB
磁盘 1 (D: E: F:) TOSHIBA MQ01ABD100 932 GB
1.1.2 操作系统
Windows10 、Ubuntu 虚拟机
1.1.3 分析软件
(含版本号)
编程工具:
R 3.5.1 Perl 5.26.1
分析工具:
ART 2.5.8 sratoolkit 2.9.4 fastqc SOAPdenovo
Quast 本地blast gffcompare Augustus
1.1.4 数据库资源
NCBI Genome:https://www.ncbi.nlm.nih.gov/genome/15?genome_assembly_id=341003
NCBI PubMed:https://www.ncbi.nlm.nih.gov/pubmed/
NCBI SRA:https://www.ncbi.nlm.nih.gov/sra/SRX257750[accn]]
UniProt:https://www.uniprot.org
EPD:https://epd.epfl.ch//index.php
1.2 研究对象
物种名称 Saccharomyces cerevisiae YJM320
Genbank Accession Number:GCA_000975885.2
1.3 方法
1.3.0 实验流程图
图 1 实验流程图
1.3.1 基因组测序模拟
分别下载YJM320的sra测序数据、gff注释数据和fna基因组数据,其中sra测序地址、gff数据地址,fna数据地址。使用art illumina软件进行序列模拟,模拟参数见表1,计算理论覆盖率。采用SRA Toolkit工具验证sra数据完整性,并从中提取fastq基因序列数据。
表 1 art illumina模拟参数设置
| l(reads的长度) | f(覆盖度) | m(片段平均长度) |
|---|---|---|
| 100 | 1 | 150 |
| 100 | 2 | 150 |
| 100 | 3 | 150 |
| 100 | 5 | 200 |
| 80 | 10 | 150 |
运行代码如下
1 | #art illumina序列模拟 |
1.3.2 序列组装
采用art illumina工具分别模拟180x和3000x两套数据,对数据进行fastqc质控分析,确认数据符合质控标准后,选择一套参考基因组数据,这里选用参考基因组数据。利用bowtie2分别进行两套数据同参考基因组的比对和samtools分析,这里我特别分析了测序深度方面的数据。利用SOAPdenovo软件对序列片段进行序列组装。利用网络平台QUAST(http://quast.bioinf.spbau.ru/)进行可视化及实际覆盖率计算。
运行代码如下
1 | #art illumina工具分别模拟180x和3000x两套数据 |
#编写R语言脚本depth.r画depth180.txt和depth3000.txt的折线图
SOAPdenovo序列组装
配置文件lib.cfg
1 | | max\_rd\_len=90[LIB]avg\_ins=180reverse\_seq=0asm\_flags=3rd\_len\_cutoff=90rank=1pair\_num\_cutoff=3map\_len=32q1=paired\_dat\_m180\_1.fqq2=paired\_dat\_m180\_2.fq[LIB]avg\_ins=3000reverse\_seq=1asm\_flags=3rank=2pair\_num\_cutoff=5map\_len=35q1=paired\_dat\_m3000\_1.fqq2=paired\_dat\_m3000\_2.fq | |
将程序挂上服务器,nohup.out作为记录工作日志
1 | nohup SOAPdenovo-63mer all -s lib.cfg -K 29 -o SOAPdenovo\_out -p 20 & |
1.3.3 同源搜索
从UniProt数据库下载所有已知真菌蛋白序列,关键词为taxonomy:fungi AND reviewed:yes,与基因序列数据fna文件进行tblastn比对,比对结果通过blast92gff3.pl程序生成注释数据gff文件。
运行代码如下
1 | #创建本地BLAST数据库 |
1.3.4 从头预测
利用augustus工具对基因序列数据fna文件进行从头预测, 提取输出文档中的有用信息进行同真菌蛋白序列进行blastp比对,从结果中获得具体的基因信息,构建基因注释数据gff文件。最后,利用gffcompare分别比对同源搜索结果、从头测序结果与原基因组
运行代码如下
1 | #创建本地BLAST数据库 |
1.3.5 启动子分析与预测
从EPD数据库中下载任意一种启动子相关的DNA元件的HMM数据,去除序列中的小写atcg和未知N,对正反双链利用bootstrap法计算pvalue和score,取99%的置信区间绘制打分图和类ROC曲线(横坐标为位点,纵坐标为该距离内的启动子/所有启动子)
1.3.5 基因组可视化
使用IGV进行可视化。基因组序列:菜单栏 Genomes→Load Genome from File =》选中fasta 格式的基因组序列文档;同源预测(blast)结果:菜单栏 Files→Load from File=》选中blast 比对结果转换来的 gff3 文档;从头预测结果:菜单栏 Files→Load from File=》选中从头预测的gff3 文档;全部加载成功后,任选一个包含注释信息的区间截图保存。
二、结果
1. 基因组测序模拟
利用art illumine HighSeq2k对不同fold的数据分别模拟测序,结果见表2
表 2 测序数据与丢失值、覆盖率的关系
| l | f | m | 丢失值 | 覆盖率 |
|---|---|---|---|---|
| 100 | 1 | 150 | 3.70E-01 | 63.0% |
| 100 | 2 | 150 | 1.37E-01 | 86.3% |
| 100 | 3 | 150 | 5.05E-02 | 95.0% |
| 100 | 5 | 200 | 7.42E-03 | 99.3% |
| 80 | 10 | 150 | 4.76E-05 | 100% |
图 2 测序覆盖度f与覆盖率的关系
分析:在fold小于10的情况下,覆盖率随fold的增大而增大,其变化幅度随fold的增大而减小,f大于10后基本可保证完全覆盖。
2. 测序组装
2.1质控分析
分别对art illumina模拟180x和3000x的数据进行fastqc质控分析
图 3 fastqc 每个碱基位点的质量评估(√) 图 4 fastqc平均序列质量(分布近似正态√)
图 5 fastqc GC含量(理论与实际相近) 图 6 fastqc 序列重复度(√)
其峰值对应GC含量与基因组相近(约38)
Fastqc结果数据符合质控标准,通过质控检验。
2.2数据比对
利用bowtie2分别进行两套数据同参考基因组的比对和samtools分析,这里我特别分析了测序深度方面的数据。结果如下
图 7 flagstat 180x
图 8 flagstat 3000x
二者的测序均在98%以上,情况符合预测。
图 9 测序深度与数量关系 上面为180,下面为3000
趋势正常(测序深度集中在45左右,符合45x),大致符合正态分布。
2.3序列组装
利用SOAPdenovo软件对序列片段进行序列组装具体结果汇总表如下:
Table 1 序列组装结果数据
| 数量 | N50 | |
|---|---|---|
| Scaffold | 265 | 258107 |
| contig | 43082(大于100bp) | 273 |
2.4覆盖率计算
利用网络平台QUAST进行可视化及实际覆盖率计算。
图 10 Quast可视化
图 11 Quast数据
利用quast分析所得的覆盖率为87.6%。
3 同源搜索
将已知真菌蛋白序列与基因序列数据fna文件进行tblastn比对,比对结果通过blast92gff3.pl生成注释数据gff文件。
4. 从头预测
利用augustus工具针对酵母基因组序列数据进行从头预测,并将结果文档同真菌蛋白序列进行blastp比对,从结果中获得具体的基因信息,补全基因注释数据gff文件。
利用gffcompare分别比对同源搜索结果、从头测序结果与原基因组,通过总结画出同源搜索、从头预测、原始基因组的对比表
表 3 同源搜索、从头预测、原始基因组的gffcompare
从头预测的结果外显子、内含子等各项情况均远优于同源建模,有理由认为从头预测优于同源建模。
5 启动子分析与预测
根据TATA box的HMM数据,对正反双链利用bootstrap法计算pvalue和score,分别绘制1号和12号染色体的打分图
图 12 染色体1全部序列的打分图
图 13 染色体12全部序列的打分图
从12号染色体的打分图中,我们也可以看出,该染色体有很长一段区域打分情况一直稳定在一个固定的值处,可以认为并没有被匹配上TATA box序列
绘制类ROC曲线(横坐标为位点,纵坐标为该距离内的启动子/所有启动子)
图 14 染色体1 +链TATAbox
图 15 染色体1 -链TATAbox
匹配得到的TATA box序列大多在基因前5000bp内
6. 基因组可视化
图 16 基因组可视化结果
从上到下依次为原基因组、同源建模结果、从头预测结果、TATA box预测结果可视化结果。